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Abstract 

.^ , A new, improved split-step backward Euler (SSBE) method is introduced and an- 

^^ I alyzed for stochastic differential delay equations(SDDEs) with generic variable delay. 

The method is proved to be convergent in mean-square sense under conditions (As- 
sumption [3lT]) that the diffusion coefficient g{x, y) is globally Lipschitz in both x and y, 
but the drift coefficient f{x, y) satisfies one-sided Lipschitz condition in x and globally 
Lipschitz in y. Further, exponential mean-square stability of the proposed method is 
investigated for SDDEs that have a negative one-sided Lipschitz constant. Our results 
show that the method has the unconditional stability property in the sense that it can 
well reproduce stability of underlying system, without any restrictions on stepsize h. 

C"^ ■ Numerical experiments and comparisons with existing methods for SDDEs illustrate 

^-v , the computational efficiency of our method. 

^ ■ AMS subject classification: 60H35,65C20,65L20. 
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- ■ 1 Introduction 

In this paper we consider the numerical integration of autonomous stochastic differential 
delay equations (SDDEs) in the Ito's sense 

dx{t) = f{x{t),x{t - T{t)))dt + g{x{t), x{t - T{t)))dw{t) (1.1) 

with initial data x{t) = ip{t),t G [— r, 0]. Here r(t) is a delay term satisfying r(t) > and 
-r := inf{t - r(t) : t > 0}, / : M'^ x M*^ — ^ M*^, ^ : M*^ x M'^ — > M'^^™. We assume that 
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the initial data is independent of the Wiener measure driving the equations and w{t) is an 
m-dimensional Wiener process defined on the complete probability space {fi,J^, {J-(}(>o,IP) 
with a filtration {J^t}t>o satisfying the usual conditions (that is, it is increasing and right 
continuous while J-q contains all P-nuU sets). 

For a given constant stepsize h > 0, we propose a split-step backward Euler (SSBE) 
method for SDDEs (11. ip as follows 

yn = yn + hf{y*^,y*J, (1.2a) 

Vn+l = y*n + 9{yl, y*n)^Wn, (1.2b) 

where Aw„ = w{tn+i) — w{tn) and for < /i < 1, 1 < g„ G Z+ 

~* ^ f ^(^n -r(t„)), t„-r(t„)<0, .^^. 

" I ^^yl-q,.+l + (1 - ^Ay*n-q^. o < t„ - r(t„) g [t„_,„, t„_g„+i). ^ • > 

For arbitrary stepsize h > 0, yn denotes the approximation of x{t) at time tn = nh,n = 
0, 1 • • ■ . We remark that fi in (11.31) depends on how memory values are handled on non-grid 
points. Generally there are two ways, the first is to use piecewise constant interpolation, 
corresponding to // = 0, and the second to use piecewise linear interpolation. In later 
development, we prefer to assume < /i < 1 to cover both cases. Also, we mention that 
the scheme (ll.2ap - (]1.2bp here is quite different from the SSBE method in [23], which will be 
explained at the end of this section. 

In (ll.2ap - (]1.2bp . y* serves as an intermediate stage value, and in order to continue the 
process, we have to solve the implicit equation (ll.2ap at every step to acquire y*. Existence 
and uniqueness of solutions to the implicit equations (ll.2ap will be discussed in section 4. 
Here, we always assume that numerical solution of (ll.2ap exists uniquely. And one can easily 
check that yn,yn is J-i„-measurable. 

The key aim in this work is to propose a new SSBE method for SDDEs with variable delay 
and its convergence and stability in mean-square sense are investigated under a non-globally 
Lipschitz condition. This situation has been investigated in [3, El El [lOl [11], [131 [El [21] for 
stochastic differential equations (SDEs) without delay. For SDEs with delay, most of previous 
work has been based on the more restrictive assumption that the coefficients /, g satisfies 
global Lipschitz and linear growth conditions, see, for example, [H [SI [151 [13, [23] • In [IS], the 
authors showed that the numerical solution produced by Euler-Maruyama (EM) method will 
converge to the true solution of the SDDEs under the local Lipschitz condition. Note that 
the proof of the convergence result in this paper is based on techniques used in [3 [18]. In [7], 
by interpreting the implicit method SSBE as the EM applied to a modified SDE the authors 
were able to get a strong convergence result. This paper, however, provides an alternative 
way to get the convergence result for SSBE. That is, by giving a direct continuous-time 
extension we accomplished the convergence proof for SSBE without considering the modified 
SDDEs. Also, in deriving moment bounds of numerical solution, due to the delay term of 
our SSBE, i.e., ^* in (ll.2ap . y* cannot be exphcitly dominated by yn as (3.25) in [7]. Starting 
with a recurrence of y* given by substituting (]1.2bp into (ll.2ap . we overcome this difficulty 
and obtained the desired moment bounds. Note that a similar approach is adopted in the 
stability analysis. 



Of course, the most important contribution of this work is to propose an improved SSBE 
method for SDDEs and to verify its excellent stability property. In [23], the authors proposed 
a SSBE method for a linear scalar SDDE with constant lag and its convergence and stability 
are studied there. It is worth emphasizing that our proposed method is a modified version 
of SSBE in [23]. The changes are in two aspects: firstly, we drop the stepsize restriction 
h = -, n E Z+ and allow for arbitrary stepsize h > 0; secondly and most importantly, 
the scheme has been modified to a new one. To see this, the two methods are applied to 
a linear scalar SDDE in section O One can observe that the second terms of /, g in the 
scheme in [23] is the numerical solution i/n-K+i (see (15. 4 p below). While the corresponding 
terms in our scheme is the intermediate stage value i/n-K. (see (15.31) below). Note that the 
modifications of the method do not raise the strong order of the numerical solution, but 
they indeed improve the stability of the method greatly. In fact, it is shown below that 
our method can well replicate exponential mean-square stability of nonlinear test problem, 
including the linear test equation as a special case, without any restrictions on stepsize h. 
The convergence and stability results of SSBE can be regarded as an extension of those in 
[3 [8] for SDEs without delay to variable delay case. This unconditional stability property 
of (ll.2ap - fll.2b|) demonstrates that the proposed method is promising and will definitely be 



effective in solving systems with stiffness in the drift term, where stability investigations are 
particularly important. 

This article is organized as follows. In next section, a general convergence result (Theorem 
12. 4p is established. In section 3, a convergence result is derived under a one-sided Lipschitz 
condition (Assumption 13 . 1 p . Section 4 and 5 are devoted to exponential mean-square stability 
property of the method. Numerical experiments are included in section 6. 

2 The general convergence results 

Throughout the paper, let | ■ | denote both the Euclidean norm in M.'^ and the trace norm(F- 
norm) in M'^^'". As the standing hypotheses, we make the following assumption. 

Assumption 2.1 The system / li.ij) has a unique solution x(t) on [— r, T] . And the functions 
f{x,y) andg{x,y) are both locally Lipschitz continuous inx andy, i.e., there exists a constant 
Lr such that 

\fix2,y2) -/(xi,?/i)p V \gix2,y2) - gix^yi)]"^ < Ln{\x2 -xip + \y2-yi\'^), (2.1) 
for all Xi,X2,yi,y2 G ^'^ with \xi\ V \x2\ V \yi\ V \y2\ < R. 

Moreover, we assume that [TSj 
Assumption 2.2 iplt) is Holder continuous in mean-square with exponent 1/2, that is 

Em)-^{s)\^<7^,\t-sl (2.2) 

and T{t) is a continuous function satisfying 

\rit) - t{s)\ < r]2\t - s\. (2.3) 



In the following convergence analysis, we find it convenient to use continuous-time approxi- 
mation solution. Hence we define continuous version y{t) as follows 

^^'' 1 yn + {t-Qf{y:,y:)+g{yl,y*n)^Mt), t e [t„, Wi),n > 0, ^ ' ' 



where Awn(t) = w{t) — w(tn)- For t G [tn,tn+i) we can write it in integral form as follows 

m ■.= yo+ [ fiy*is),y*{s))ds+ [ giy*is),y*{s))dw,, (2.5) 

Jo Jo 

where 



y*is) := Y^ l{i„<s<t„+i}2/^, y*{s) ■= Y^ l{t„<s<t„+^}yn- (2-6) 

n=0 n=0 

It is not hard to verify that y{tn) = yn, that is, y(t) coincides with the discrete solutions at 
the grid-points. 

In additional to the above two assumptions, we will need another one. 

Assumption 2.3 The exact solution x{t) and its continuous-time approximation solution 
y{t) have p-th moment bounds, that is, there exist constants p > 2,A > such that 



E 



sup |a:(t)|^ 

0<t<T 



VE 



sup imr 

0<t<T 



<A. 



(2.7) 



Now we state our convergence theorem here and give a sequence of lemmas that lead to 
a proof. 

Theorem 2.4 Under Assumptions \2. 1^2.2^2. 3\ if the implicit equation lil.2a\) admits a unique 



solution, then the continuous-time approximate solution y{t) (2^ will converge to the true 
solution of ( (i.ij) in the mean-square sense, i.e., 

E sup |y(t) -x(t)|^ ^0, as h-^0. 

0<t<T 

We need several lemmas to complete the proof of Theorem 12. 4[ 
First, we will define three stopping times 

p^ = inf{t>0:|x(t)|>i?}, r^ = inf{t>0: \y{t)\ > R ,oi \y*{t)\ > R}, aR = pnATR, 

where as usual inf is set as oo (0 denotes the empty set). 



Lemma 2.5 Under Assumption \2.1\ \2.2[ there exist constants Ci{R), C2{R) such that for 
s G [t„,t„+i) and h < 1 



El{.<.«}|y(s)-y*(s)|2<Ci(/2)/i, 

^ks<.Jy{s - t{s)) - r(s)P < C2{R)h. 



(2.8) 
(2.9) 



Proof. For s G [t„,t„+i), by definition oi y{s) and y*{s), 

y{s)-y*{s) = yn + {s-tn)f{y*n,y*n)+9{yn,yn)^Wn{s)-y*n 

= {s-K+{)f{yl,rn)+9{ylyn)^Wn{s)- (2.10) 

Noticing that for |x| V |'y| < R 

|/(x,y)r < 2|/(x,y)-/(0,0)p + 2|/(0,0)p 

< Kr(1 + |xP + ||/P), (2.11) 

with Kr = 2 maxjLij, |/(0, 0)|}. Using hnear growth condition of g and moment bounds in 
(13. 7p . we have appropriate constant Ci{R) so that 

El{.<.«}|y(s) - y*is)\' < 2Knh\l + ny*^ + El^H + 2i^/i(l + ny*S + ^1^1') 

< Ci{R)h. 

As for estimate (12.91) . there are four cases as to the location of t„ — r(t„,) and s — t{s): 

• 1) tn-T{tn) < 0,S-r(s) < 0, 

• 2)t„-r(t„)>0,s-r(s)>0, 
•3) t„-r(t„) <0,s-r(s) > 0, 

• 4) t„-r(t„) >0,s-r(s) < 0. 

Noticing that the delay t{s) satisfies Lipschitz condition (12. 3p . one sees that 

\S - t{s) -tn + T{tn)\ < {j]2 + l)h. (2.12) 

In the case 1), combining Holder continuity of initial data (12. 2p and (I2.12p gives the desired 
assertion. In the case 2), without loss of generality, we assume s— r(s) G [tj, tj+i), t„ — r(t„) = 
(1 — ij)tj + /itj+i G [tj,tj+i), "i > j > 0. Thus we have from (ll.2ap and (13.80 that 

y(s - t{s)) - y*{s) = y^ + {s- r{s) - U)f{y*, y*) + g{y*, y*)Awi{s - t{s)) 

-{'^ - f^)y*j - f^yj+i 
= {s- t{s) - U+i)fiyh y*) + giy*, y*)Aw,{s - r(s)) 

+(1 - f^)iyi - y*j) + ^(^r - y*j+i) 

= {s- t{s) - U+i)f{yl y*) + 9{yl y*) Ai/;.(s - t{s)) 

i-l 

+(1 - /") Zl [hf(yl+i^yl+i) + g{ylyl)AM 

k=j 
i-l 

+/" E [hf{yl+i,yl+i) + 9{yl,yl)^w,] , (2.13) 

where as usual we define the second summation equals zero when i = j + 1. Noticing from 
(I2.12P that i—j < r/2 + 1? and combining local linear growth bound (12.111) for /, global linear 
growth condition for g and moment bounds (13.71) . we can derive from (I2.13P that 

El{,<^^}|r(s) - y{s - r(s))|2 < C2{R)h. 



In the case 3) and 4), using an elementary inequality gives 



<2El{,<,^}|r(s) - Vm' + 2El{,<,,}|y(0) - vis - r(s))|2. 

Then combining this with results obtained in case 1) and 2) gives the required result, with 
6*2 (-R) a universal constant independent of h. 



Lemma 2.6 Under Assu'mption \2.1\ \2.2\ for stepsize h < 1, there exists a constant Cr such 
that 

E[ sup \y{t AaR)-x{tA a^)]^] < CrH, 

0<t<T 

with Cr dependent on R, hut independent of h. 
Proof. For simplicity, denote 

e{t):=y{t)-x{t). 
From (11. ip and fl2.5p . we have 



E 
E 



sup \y[s A aR) - x[s A (Tr) 

0<s<t 



sup \e[s A o-rW = E 

0<s<t 

sup / f{y*{r),y*{r))- f{x{r),x{r-T{r)))dr 

0<s<t Jo 
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iAan 



9iy*ir),y*{r)) - g{x{r),x{r - r(r)))dw(r) 







< 2TE/ \f{y*{s),y*{s))-f{x{s),x{s-T{s)))fds 

Jo 

fsAan, 



-2E 



sup 

0<s<t 



9{y*{r),y*{r)) - g{x{r),x{r - T{r)))dw{r) 



rtAa-R 

< 2(r + A)LrE / \y*{s) - x(s)p + iris) - xis - r(s))pds, (2.14) 

Jo 

where Holder's inequality and the Burkholder-Davis-Gundy inequality were used again. Us- 
ing the elementary inequality \a + 6p < 2|ap + 2|6p, one computes from (I2.14p that 

E[sup \eisAaR)\^] 

0<s<t 

ptAa-R 

< 4{T + 4)LrE \y*(s)-y{s)\^ + \y{s)-x{s)\Ms 

Jo 

rtAo-R 

+4(T + 4)LrE / \r{s) - y{s - t{s))\^ + \y{s - t{s)) - x{s - T{s))\'ds 
Jo 



< 8(T + 4)Lr [ E[ sup \yir A a^) - a;(r A a^)|']ds 

Jo 0<r<s 

+4(T + 4)L«E / |y*(s)-y(s)|2ds 
Jo 



rtAan 

+4{T + 4)LrE / |^*(s)_^(s_r(s))pds, (2.15) 

where the fact was used that \y{s — r(s)) — x{s — r(s))p < sup \y{r)) — x(r)p. By taking 

0<r<s 

Lemma 12751 into account, we derive from fl2.15p that, with suitable constants Cr,Cr 
E[sup \e{sAaR)\^] < 8(T + 4)Lr / E[ sup \y{r A cr) - x{r A(TR)\^]ds 

0<s<t Jo 0<r<s 

+4(T + A)TLRCi{R)h + 4(T + A)TLRC2{R)h 
= Cr I E[snp \e{r A aR)\^]ds + CRh. (2.16) 

Jo 0<r-<s 

Hence continuous Gronwall inequahty gives the assertion. 

Proof of Theorem \2.4\ Armed with Lemma 12.61 and Assumption 12. 3[ the result may be 
proved using a similar approach to that in [TJ Theorem 2.2] and [HI Theorem 2.1], where 
under the local Lipschitz condition they showed the strong convergence of the EM method 
for the SODEs and SDDEs, respectively. 



Remark 2.7 Under the global Lipschitz condition and linear growth condition (cf fTT^ ), we 
can choose uniform constants Ci{R), C2{R), Cr in previous Lemma \27^2. 6\ to be independent 
of R. Accordingly we can recover the strong order of 1/2 by deriving 

E[ sup \y{t) -x{t)\^] < Ch, 

0<t<T 

where C is independent of R and h. 

3 Convergence with a one-sided Lipschitz condition 

In this section, we will give some sufficient conditions on equations (11. ip to promise a unique 
global solution of SDDEs and a well-defined solution of the SSBE method. We make the 
following assumptions on the SDDEs. 

Assumption 3.1 The functions f{x,y) are continuously differentiable in both x andy, and 
there exist constants 71,72,73,74, such thatWx,y, Xi,X2, yi,y2 G ^'^ 

{x2-xij{x2,y)- fixi,y)) < -fi\x2 - xi\'^ , (3.1) 

\f{x,y2)- f{x,yi)\ < 72I1/2-Z/1I, (3.2) 

\g{x2,y2)- gixi,yi)\'^ < 73|a;2 - a;i|^ + 74I2/2 - yiP- (3.3) 



The inequalities (13. ip .( 1312]) indicate that the first argument x of / satisfies one-sided Lips- 
chitz condition and the second satisfies global Lipschitz condition. It is worth noticing that 
conditions of the same type as (13.1 p and (13.20 have been exploited successfully in the anal- 
ysis of numerical methods for deterministic delay differential equations (DDEs) (see [3] and 
references therein). As for SDEs without delay, the conditions (13. ip and (13. 3 p has been used 
in [3 El Ilia [21]. 

We compute from (E1])-([S3D that 



{xj{x,y)) = {xj{x,y)-f{0,y)) + {xj{0,y)-fi0,0)) + {xj{0,0)) 
< (7i + l)N^ + ^72|z/P + ^|/(0,0)p, 

\g{x, y)\^ < 2\g{x, y) - g{0, 0)\' + 2\g{0, 0)p < 2^,\x\^ + 2^^' + 2\g{0, 0)\'. 
On choosing the constant K as 

ir = max |7i + 1,273, ^72, 274, ^1/(0, 0)|2, 21^(0,0)1' 

the following condition holds 



(3.4) 
(3.5) 



X 



^f{x,y) V \g{x,y)\^ < K{1 + \x\^ + \y\^), \/x,y G 



(3.6) 



In what follows we always assume that for Vp > the initial data satisfies 



E 



E sup \ij{s] 

-T<S<0 



< oo. 



Theorem 3.2 Assume that Assumption \3.1\ is fulfilled. Then there exists a unique global 
solution x{t) to system U.l\) . Morever, for any p > 2, there exists constant C = C{p,T) 



E 



sup |a;(t)|^ 
o<t<r 



<cii+Emn. 



Proof. See the Appendix. 



Lemma 3.3 Assume that f,g satisfy the condition Ii3.6\) and h < 1 is sufficiently small , 
then for p > 2 the following moment bounds hold 



E 



sup \yn 

0<nh<T 



2p 



VE 



sup \r{t)\'p 

0<t<T 



VE 



sup \y{t)\ 

0<t<T 



2p 



< A. 



(3.7) 



Proof. Inserting fll.2bp into (ll.2al) gives 

y* = 2/*_i + hf{y*^, y*J + g{yl_i, yl_i)Awn-i, n > 1. 
Hence 

IVn - hf{y*n,y*n)\^ = K^i + ^(^^.i, C_i) Aw„,_i |^ 



(3.^ 



Expanding it and employing (13. 6p yields 



,*|2 



2Kh{l + \y*J+\y*S) < \y:_,\' + 2{y:_„ giy:_,,y*n-i)^w^~i) 

+ \9{yn-l,yn-l)^Wn-l\'^ . 



(3.9) 



By definition of y*, one obtains |y*p < |?/*p+maxo<i<n-i |y*P + ||'?/'p. Taking this inequality 
into consideration and letting h < ho < 1/(47^"), we have from (13. 9 p that 



(1 - AKh)\y:\' < \y:_,\' + 2Kh{l + max {y*]' + ||^|n 

0<«<n— 1 

Denoting a = 1/(1 — AKho), one computes that 

|Z/;^P < |y;„i|^ + 6fs:a/i max |y*|2 + 2fs:a/i + 2Js:a/i||^||2 

0<i<n—l 



(3.10) 



n-1 



(3.11) 



By recursive calculation, we obtain 



n-l 



\yl? < \yl? + QKah'S^ m&x\y*? + 2KaT + 2KaT 



i=0 



0<i<i 



n— 1 n— 1 

-2a ^ (y*, (7(2/*, y;)A«;,) + a ^ |(7(2/;, y*)A 

j=0 j=0 



w^ 



Raising both sides to the power p gives 

|^*|2p ^ 5P-1 I |^*|2p ^ (6i^a/i)PnP-i ^ max \y*\^P + [2KaT + 2KaT\\ijfY 



n-l 



j=0 



0<i<j 



-(2a) 



n-l 



.i=o 



a^n^^ 



n-l 
j=0 



Thus 



E max \y* 

l<n<M '^" 



* |2p 



M-1 



< 5^-^ <^ E|y*r^ + {6KafTP-^hE ^ max li/*!^^ + E (2i^aT + 2i^aT||^|p) 



i=o 



+ (2a)^E max 

l<n<M 



n-l 



^{yj,9iy*j.yj)^Wj) 

.3=0 



M-1 



+ a^M^-'E Y, \g{y*,y*)Aw,\'p \ . (3.12) 

i=o 



Here 1 < M < A^, where A^ is the largest integer number such that Nh < T. Now, using 
the Burkholder-Davis-Gundy inequality (Theorem 1.7.3 in [17]) gives 



E max 

l<n<M 



n-l 


p 


'M-1 


y](y*,9ijj;,y*)Aw,) 


<CpE 


"Yly^'giyhm'h 


j=o J 




. i=o 



p/2 



< Cp{Khy/^MP/^~^E 



■M-l 



*|2 I |~*|2\P/2 



Ei^;i'(i+i^;i'+i^;n 



- 2 ^ 



j=0 

M-l 



L i=o 



(3.13) 



Noticing that 



E\v*A'^P < E max \v*\^p + EII^II^p 



0<i<j 



(3.14) 



inserting it into (I3.13p . we can find out appropriate constants C = C(j),K,T) such that 



E max 

0<n<Af 
A/-1 



n-l 



J2{yh9iyhyj)Aw,) 



j=0 



< C/iVEm£ 
■^^^ 0<i 



i=o 



max \yi 

<3 



\2p 



C(E 



2p 



1). 



(3.15) 



At the same time, noting the fact y^^iin ^ ^tn and Aw^ is independent of J-'i^, one can 
compute that, with C = C{p,T) a constant that may change hne by hne 



Af-l 



M-l 



^Y.\9iyhy*,)^^3\'' < J2^\9{yhy*)\''n^w, 



|2p 



i=o 



i=o 



A/-1 



j=0 

Af-l 

< C'/iP-^(E||^fP + l) + C'/iPy Emax|y*pP. (3.16) 



i=o 



By definition (ll.2al) . one sees that 



\y*o-hf{y*o,yo)\' = \yo\'. 

Then using a similar approach used before, we can find out a constant cq = co{p,K) to 
ensure that 

E||/*pP<Co(E||^||2p + l)<oo. (3.17) 

Inserting (I3.15p . (l3.16|) into (13.121) and considering (13.171) and /i < 1 , we have, with suitable 
constants C = C'{p, K, T), C" = C"{p, K, T) 

E max ||/:PP < E|y*pP + E max |y:P^ 



A//-1 



< C'iEUf^ + 1) + C'/i X E max. ly^'^ 



i=0 



o<j<i 



(3.18) 
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Thus using the discrete- type Gronwall inequahty, we derive from fl3.18p that E [supo<„/j<T Iz/j^P^] 
is bounded by a constant independent of N. Then by considering the elementary inequahty 
\fix + (1 — /i)ypP < /i|xp^ + (1 — fi)\y\'^P, boundedness of E [suPq<^<2" IJ'*(^)P^] is immediate. 
To bound E [suPq<^<^ |y(t)p^], we shah first bound E [suPq<„^<7. ll/nP^]- From f ll.2bp . 
we have 



E 



sup \yr 

0<nh<T 



|2p 



< 22P-^<'E 

< 22p-i I E 



sup \yj 

0<nh<T 



+ E 



sup Igivn, y*n)^w, 

0<nh<T 



\2p 



sup |2/„ 

0<nh<T 



2p 



N 



j=0 



Now (13.161) and bound of E [supo<„,/i<T |y^P^] gives the bound of E [supo<„/i<r Iz/nP^] • 
To bound E [suPq<«2. |y(t)p^]^ we denote by rit the integer for which t G [tnt,tnt+i] 
definitions of flL^ and (^, for t > 0, 

y{t) = ynt + {t-tn,)f{y*r,^,ynj+9{ynt,y*nj^^nt{t) 
= Vnt + l{y*n, - Vnt) + giVnt^ Vl, ) ^Wn^ (t) 

= (1 - i)ynt + ivL + givL^ynJ^WnM, 



By 



(3.19) 



where 7 = (t — t„,J//i < 1. Thus 



E 



sup \y{t)\'^ 

0<t<T 



sup |2/„ 

0<nh<T 



2p 



sup \giyl^,ylJAwn,it)\ 

0<t<T 



1 -7)E 

2p 



< 22P"^<J7E 
+E 
Using Doob's martingale inequality [TTJ Theorem 1.3.8], we derive that 



sup \yn 

0<nh<T 



2p 



(3.20) 



E 



sup \giy*y*)AwnM\ 



2p 



0<t<T 



N 



< 



Ee 



n=0 



sup \giyl,y^)Awnis] 

0<s<h 



|2p 



< 



2p 



2p 



■s2pN 

- J2^[\g{y:,y:)^wM\"'']. (3.21) 

^ n=0 



Thus the last term in (13.201) is bounded by considering (I3.16P and bounds of E[supQ<,„/^<y |y, 
E[supo<„/j<T l^nP^]- Now boundedness of E[suPq<^<7. Ii/(^)P^] follows immediately. 



* |2pl 
n\ V 



Lemma 3.4 Under Assu'mption \3.1l ^7(71 + 72)^ < 1? the implicit equation in U.2a\) admits 
a unique solution. 



Proof. Let /(c) := f{c,iJ,c+ (1 — jj,)b), then the implicit equation (ll.2ap takes the form as 

c = hf{c) + d = hf{c, /ic + (1 — fi)b) + d, 
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where at each step, < /i < l,b,d are known. Observing that 

(Cl-C2,/(Ci) -/(Ca)) = (Ci -C2,/(Ci,/iCi + (l-/i)6) -/(c2,/iCi + (1 -;U)6)) 

+ (Cl - C2, /(C2, /iCi + (1 - /i)&) - /(C2, /iC2 + (1 - /i)&)) 

< 7l|Cl -C2p + /i72|Ci -C2p 

< (71 +72)|Cl -C2|^ 

the assertion follows immediately from Theorem 14.2 of [6]. 

Corollary 3.5 Under Assumption \2.0i3.1[ if (71 + 72)/^ < 1, then the numerical solution 
produced by i[1.2a\) - D.2b\) is well-defined and will converge to the true solution in the mean- 
square sense, i.e., 

E sup |y(t) - x(t)|^ -)■ 0, as /i -)■ 0. 
o<t<r 

Proof. Noticing that Assumption 13.11 implies Assumptions I2.1|2.3I by Theorem 13.21 and 
Lemma 13.31 ^"^^ taking Lemma 13.41 into consideration, the result follows directly from The- 
orem |2]H 

Remark 3.6 We remark that the problem class satisfying condition \2.3\) includes plenty of 
important models. In particular, stochastic pantograph differential equations (see, e.g., J^) 
with T{t) = (1 — q)t, < g < 1 and SDDEs with constant lag fall into this class and therefore 
corresponding convergence results follow immediately. 

4 Mean-square stability with bounded delay 

In this section, we will investigate how SSBE shares exponential mean-square stability of 
general nonlinear systems. In deterministic case, nonlinear stability analysis of numerical 
methods are carried on under a one-sided Lipschitz condition. This phenomenon has been 
well studied in the deterministic case ([21 Ej and references therein) and stochastic case 
without delay [3 [3 [9l [121 l2l]- In what follows, we choose the test problem satisfying 
conditions ( I3.ip -f l33|) . Moreover, we assume that variable delay is bounded , that is, there 
exists r > 0, for 1 < K G Z+, < 5 < 1 

0<r(t)<r, T = {K-5)h. (4.1) 

We remark that this assumption does not impose additional restrictions on the stepsize h 
and admits arbitrary large h on choosing n = 1 and < 5 < 1 close to 1. To begin with, 
we shall first give a sufficient condition for exponential mean-square stability of analytical 
solution to underlying problem. 

Theorem 4.1 Under the conditions Ii3. 1\) . ^37^) . / 1 5*. 3\) and (^TTj), and with 71,72,73,74 obey- 
ing 

/3:=27i + 272 + 73 + 74<0, (4.2) 

any two solutions x{t] ip) and y{t; (p) with E||^p < 00 and E||</)p < 00 satisfy 

E\x{t) - i/(t)p < E||0 - ipfexp{-u+t}, 
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where u^ G (0, — /3] is the zero of C{u) = v + Pi+ (32 exp{z/r}, with Pi = 271 + 72 + 73, /32 = 

72 + 74- 

Proof. By Ito formula, we have 

E\x{t + 6)- y{t + 5)|2 - E\x{t) - y(t)|2 
t+s 

2E{xis) - vis), f{x{s), x{s - ris))) - /(y(s), y{s - T{s))))ds 
t 

rt+5 

+ / E\g{x{s),x{s - t{s))) - 9{y{s),y{s - r(s)))|Ms 

/t+5 r-t+5 

E\x{s) - y{s)\^ds + 74 / E\x{s - t{s)) - y{s - T{s))\^ds 

i-t+S 

+2 j E{x{s) - y{s), f{y{s),x{s - t{s))) - f{y{s), y{s - T{s))))ds 

rt+S j-t+S 

< Pi E\x{s) - y{s)\'^ds + P2 sup E\x{r) - y{r)\'^ds. (4.3) 

Jt Jt re[s-T,s] 

Letting u(t) = E\x(t) — y(t)p and noticing that u{t) exists for t > — r and is continuous, we 
derive from (14 .Sp that 

D+M(t) < Piu{t) + P2 sup u{s), 

S€[t-T,t] 

where the upper Dini derivative D^u(t) is defined as 

_^ ^^^ ,. u{t + 6) - u{t) 



D u(t) := hmsup 

Using Theorem 7 in [2] leads to the desired result. 

Based on this stability result, we are going to investigate stability of the numerical 
method. 

Theorem 4.2 Under the conditions ( TO) . fS^) . ( TOI) and ^J\ ), if P <0, then for allh > 0, 
any two solutions Xn,Yn produced by SSBE U.2a\) - n.2b\) with EW^pW^ < 00 and IE||0|p < 00 
satisfy 

E|X„ - y„|^ < E||0- V^|pexp{-z/^ra/i}, as n ^ 00, 

where u'^ > is defined as 

'^ 2{K+l)h \1 + h-f2 + hys + h-fj ^ ' 

Proof. Under /3 < 0, the first part is an immediate result from Lemma 13.41 For the 
second part, in order to state conveniently, we introduce some notations 

w: = x:- y;. A/: = f^x^K) - f{y:.Y:). ^9: = giK.K) - giy^X). (4.5) 

From (13. Sp . we have 

W: = W:_i + HAf: + Ag:_iAwn^i. (4.6) 
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Thus 

Taking expectation and using (13.31) yields 

E\w:\' - 2hE{w:,Af:) < (i + hr,)E\w:^,\' + h^,nx:_, - y:_,\'. (4.?) 

Now using the Cauchy-Schwarz inequahty and conditions f l3.ip -( l3l2l) . we have 

2E{w:,Af:) = 2E{w:jix:,x:)-fiY:,x:)) 
+2E{w:j{Y:,x:)-f{Y:,Y:)) 

< 2^iE\W: 1 2 + 272EI PV; \\K-K\ 

< (27i + 72)E|iy:|2 + 72E|l:-F;r. 

Inserting it into (14. 7p gives 

(1 _ 2/,7i - h^2)E\x: - y:\' < (1 + h^s)E\x:_, - y:_,\' 

+h-f,E\x:_, - ?:_,{' + h^,E\x: - y:\\ (4.8) 

Here we have to consider which approach is chosen to treat memory values on non-grid 
points, piecewise constant interpolation (/x = 0) or piecewise linear interpolation. In the 
latter case, let us consider two possible cases: 
• If r{tn) = jlh, < /i < 1, then 



(4.9) 



E\x: - f;p = E\~^xu + (1 - fi)x: - ~^y:_, - (i - ~^)y:\' 
< mxu - y^i? + (1 - m\K - y:\\ 

Inserting (14. 9p . we derive from (14. 8 P that 

[l_2/i7i-(2-/i)/i72]E|X:-F;p 

< (1 + /173 + M72)E|x:_i - f;_i|2 ^ /,^^E|x:_i - y:_,\\ 

Hence using the fact /3 < in (14. 2 p gives 

^l^'^-^-l ^ l-2/.7i-(2-/I)/.72 ----^l"<»-i'^'''-"'^' 

< ^ + ^^' ^ ^^' t ^^' max E|X;-F;r. (4-10) 

1 — 2/171 — /l72 n-K-l<i<n-l 

• If r(t„) > /i, it follows from (|MD and /3 < that 

JM. A„ — Y„\ < ; ; max E\X, — K- . 4.11 

1 — 2/?-7l — /l72 7i-K-l<'t<n-l 

Therefore, it is always true that inequality (14. lip holds for piecewise linear interpolation 
case. Obviously (14. lip also stands in piecewise constant interpolation case. 
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Further, from fll.2ap one sees 



\X* - Y* - hifiX*,X*) - fiY*,Y*))\' = \Xo - Fop. 
Using a similar approach as before, one can derive 

Denote 

A. := '^ "'' + ''"%+"" . (4.13) 

1 — 2^71 — A172 

Noticing that (3 < 0, one can readily derive < /3/i < 1, we can deduce from (14.111) and 
gl2D that 

I ^~2 I I 1 71 — 2 

Here \_x\ denotes the greatest integer less than or equal to x. 

Finally from ([L2b]), we have for large n such that 23±2i + ^^^ ln/3,, < ^(^^ 

E|x„-r„p < ii + h^s)E\x:_,-Y:^,\' + h^,E\x:_,-Y:^,\' 

n~2 n — K, — 2 

< (1 + /i73)/3r'^ii^ - 0f + /^74/37^i?ii^ - 0ir 

n-K-2 

< E||0-V^f exp{-z/+n/i}, (4.14) 

where z/^ is defined as in (14.41) . 

The stability result indicates that the method (ll.2al) - fll.2bp can well reproduce long-time 



stability of the continuous system satisfying conditions stated in Theorem 14.11 Note that the 
exponential mean-square stability under non-global Lipschitz conditions has been studied in 
[8] in the case of nonlinear SDKs without delay. The preceding results can be regarded as 
an extension of those in [8] to delay case. 

5 Mean-square linear stability 

Although the main focus of this work is on nonlinear SDDEs, in this section we show that 
the SSBE (ll.2al) - (11.2bp has a very desirable linear stability property. Hence, we consider the 
scalar, linear test equation [151 ES] given by 

dx{t) = {ax{t) + hx{t - r))dt + (cx(t) + dx{t - T))dw{t). (5.1) 

Note that (15. ip is a special case of (II. ip with r(t) = r, and satisfies conditions (I3.ip - (l3.3p 
with 

7i = a, 72 = l&l, 73 = c^ + |crf|, 74 = rf^ + |cd|. 
By Theorem 14. 11 (15. ID is mean-square stable if 

a<-\h\-\{\c\ + \d\f. (5.2) 
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For constraint stepsize h = t/k, 1 < k G Z'^, i.e., S = in (14.11) . the SSBE proposed in our 
work applied to (15.1 p produces 

Vn+l =yn+ [Cy*n + %:_«]Aw„. 

In [23], the authors constructed a different SSBE for the hnear test equation (15.11) and their 
method apphed to (15.11) reads 

Zn = Zn + h[aZ^ + hZn-K+l\, ,^^\ 

Zn+1 = Zl+ [CZI + dZn-K+l]^Wn. 

The stabihty results there [231 Theorem 4.1] indicate that under (15. 2p the method (15.40 can 
only preserve mean-square stability of (15. ID with stepsize restrictions, but the new scheme 
( 15. 3p exhibits a better stability property. 

Corollary 5.1 For the linear equation \5. 1\) . if li5.S\) holds, then the SSBE Ii5.3\) is mean- 
square stable for any stepsize h = r/n, 1 < k G Z"*". 

Proof. The assertion readily follows from Theorem 14. 2[ 

Apparently, the SSBE (15.30 achieves an advantage over (15. 4p in stability property that 
the SSBE (15. 3p is able to inherit stability of (15. ip for any stepsize h = t/k, 1 < k G Z+. If 
one drops the stepsize restriction h = -,k E Z+ and allow for arbitrary stepsize h > 0, one 
can arrive at a sharper stability result from Theorem 14.21 



IS 



Corollary 5.2 For the linear equation / TO]) , if [5^) holds, then the SSBE M.2a\) - lil.2h\) 
mean-square stable for any stepsize h > 0. 

6 Numerical experiments 

In this section we give several numerical examples to illustrate intuitively the strong conver- 
gence and the mean-square stability obtained in previous sections. 

6.1 A linear example 

The first test equation is a linear Ito SDDE 

dx{t) = {ax{t) + bx{t - l))dt + (cx(t) + dx{t - l))dw{t), 

x{t) = 0.5, tG[-l,0]. ^ ' 

Denoting y^ as the numerical approximation to x^^'{tN) at end point t]\[ in the i-th simulation 
of all M simulations, we approximate means of absolute errors e as 



^-^Y:\y'S-y'"M 



j=l 
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Figure 1: loge with tjv = 1 versus logh for Example I (left) and Example II (right). 



Table 1: Numerical results for Example II and III with t^ = 8. 







Example II 






Example III 




h 


EM 


SSBE (lOl) 


SSBE (lO)) 


EM 


SSBE (lOl) 


SSBE (lOl) 


2-^ 
2-6 
2-5 

2-4 

2-3 


0.0008 
0.0013 
0.0021 
0.0034 
0.0086 


0.0011 
0.0016 
0.0029 
0.0058 
0.0148 


0.0008 
0.0013 
0.0019 
0.0027 
0.0038 


0.0014 
0.0025 
0.0058 
0.2744 
6.1598e+010 


0.0020 
0.0036 
0.0070 
0.0157 
0.0628 


0.0014 
0.0023 
0.0035 
0.0053 
0.0078 

















In our experiments, we use the SSBE (15. 3p to compute an "exact solution" with small stepsize 
h = 2"^^ and M = 5000. We choose two sets of parameters as follows 

• Example I: a = —2, b = l,c = d = 0.5; 

• Example II: a = —6, b = 3,c = d = 1. 

• Example III: a = -20, b = 12,c = 2,d = 1. 

In Figure [H computational errors e versus stepsize h on a. log-log scale are plotted and 
dashed lines of slope one half are added. One can clearly see that SSBE (15. 3p for linear 
test equation (16.11) is convergent and has strong order of 1/2. In Table [H computational 
errors e with tjsf = 8 are presented for the well-known Euler-Maruyama method [18], the 
SSBE method (15.41) and the improved SSBE method (15. 3p in this paper. There one can find 
that the improved SSBE method (15.31) has the best accuracy among the three methods. In 
particular, for Example III with stiffness in drift term (i.e., a = —20), when the moderate 
stepsize h = 1/8 was used, the Euler-Maruyama method becomes unstable and the two SSBE 
methods still remain stable, but with the improved SSBE (15.31) producing better result. 

To compare stability property of the improved SSBE and SSBE in |23], simulations by 
SSBE (15.31) and (15.41) are both depicted in Figure [2|, [31 There solutions produced by (15.31) 
and (15.41) are plotted in solid line and dashed line, respectively. As is shown in the figures, 
methods (15.31) and (15. 4p exhibit different stability behavior. One can observe from Figure 
[2] that (15.30 for Example II is mean-square stable for h = 1, 1/2, 1/3, 1/4. But (15.40 is 
unstable for h = 1, 1/2. For Example III, the improved SSBE (15. 3p is always stable for 
h = 1, 1/4, 1/6, 1/10, but ( 15. 4p becomes stable when the stepsize h decreases to h = 1/10. 
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The numerical results demonstrate that the scheme (15. 3p has a greater advantage in mean- 
square stability than f l5.4p . 




20 25 30 




10 20 30 40 50 




Figure 2: Simulations for (16.1 p with a = —Q,b = 3,c = d = 1. Upper left: /i = 1, upper 
right: h = 1/2, lower left: h = 1/3, lower right: h = 1/4. 



6.2 A nonlinear example 

Consider a nonlinear SDDE with a time-varying delay as follows 

dx{t) = [-Ax{t) - 3x^{t) + x{t - r(t))] dt + [x{t) + x{t - r(t))] dw{t),t > 0, 

x{t) = i, te[-i,o], 



(6.2) 



where r(t) = j^. Obviously, equation (16. 2p satisfies conditions (I3.ip -( l3l3|) in Assumption 
13. H with 7i = —4,72 = 1,73 = 74 = 2. Thus 271 + 272 + 73 + 74 = — 2 < and the 
problem is exponentially mean-square stable. As is shown in Figure HJ the SSBE (15. 3 p can 
well reproduce stability for quite large stepsize h = 1,2, 5. This is consistent with our result 
established in Theorem 14.21 



Appendix 



Proof of Theorem \3.S\ Since both / and g are locally Lipschitz continuous. Theorem 3.2.2 
of [IE] shows that there is a unique maximal local solution x{t) on t G [[0,poo[[, where the 
stopping time p^ = inf{t > : \x{t)\ > R}. By Ito's formula we obtain that for t > 



tApR 

2 / x{sff{x{s),x{s-T{s)))ds 





\xit A pn)\' = mO)\ 

rt^pR 

+2/ x{s)^ g{x{s),x{s - T{s)))dws + 
Jo JO 



\g{x{s),x{s-T{s)))\ ds 




O 2 




Figure 3: Simulations for (16. ip with a = —20, b = 12,c = 2,d = 1. Upper left: h = 1, upper 
right: h = 1/4, lower left: h = 1/6, lower right: h = 1/10. 




Figure 4: Simulations for (16. 2p by SSBE (15. 3 p using various stepsizes. 
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< \tp{0)\^ + 3K (1 + |x(s)|2 + \x{s - r{s))\')ds 

Jo 



ftApR 

+2 x{s)^g{x{s),x{s-T{s)))dWs 

Jo 

where the condition fl3.6D was used. Thus 



(6.3) 



sup|x(sAp^)P < mO)\^ + 3K [ {1 + 2 snp \xirApn)\^ + m^)ds 

0<s<t Jo 0<r<s 



SAPR 



+2 sup / x(r) g{x{r),x{r — T{r)))dwr. 

0<s<tJo 

Now, raising both sides of (16. 4p to the power p/2 and using Holder's inequahty yield 

sup |x(sApH)r< 3^/2-1 {|^(0)|P 

0<s<i 

+ (3ir)P/2(3T)P/2-i /"(I + 2^/2 g^p \a;(r A pr)\p + \\^l:\\P)ds 

Jo 0<r<s 



(6.4) 



.2P/2 



sup 

0<s« 



sApR 



x{r)'^g{x{r),x{r — T{r)))dwr 



p/2- 



(6.5) 



By the Burkholder-Davis-Gundy inequality [T7], one computes that, with ci = Ci(p, T), 



E 



sup \x{sApr)\p 

0<s<t 



< ci <; 1 + E||^f + / E sup \x{r A pR)fds 

Jo 0<r<s 

tApR -\ P/4 ■ 



+E 



\x{s)\ \g{x{s),x{s — t{s)))\ ds 



(6.6) 



Next, by an elementary inequality, 

-tApR 



E 



< E 







x{s)\ \g{x{s),x{s-T{s)))\ ds 

■tApR 

ip |x(sApij)p / 

_0<s<t Jo 



p/4 



< — E 

- 2ci 

1 ^ 

< — E 

- 2ci 



sup \x{s ApRJf 

0<s<t 

sup \x{sApr)\p 

0<s<t 



sup |x(sApij)p / \g{x{s),x{s - T{s)))\'^ds 

0<s<t Jo 

„ ptApR 

+ Atp/^-^EJ \gixis),xis-Tis)Wds 

pt 
^ ^(3y^p/2-ij^p/2 / (1 + E sup |x(rAp^j|P + E||^f)ds. 

2 Jo 0<r<s 

Inserting it into (16. 6p . for proper constants C2, C3 we have that 

E sup |x(s Apr)!*' < C2(l + E||^f) +C3 / E sup \x{r ApR)Y'ds. 

0<s<t Jo 0<r<s 
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The Gronwall inequality gives 

E sup \x{sApR)\P<C2{l + E\\i/j\\P)e'''^. (6.7) 

0<s<T 

This iinphes 

i^^PW < T} < C2(l + E||^r)e^^^. 

Letting /2 — )• oo leads to 

lim F{pR <T} = 0. 

Since T > is arbitrary, we must have p/j — )■ oo a.s. and hence poo = oo a.s. The existence 
and uniqueness of the global solution is justified. Finally, the desired moment bound follows 
from (16. 7p by letting i? — )■ oo and setting C = c^e^''"'^ ■ 
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